NASA 

Contractor  Report  195476 


Army  Research  Laboratory 
Contractor  Report  ARL-CR— 237 


An  Efficient  Numerical  Procedure  for 
Thermodydrodynamic  Analysis 
of  Cavitating  Bearings 


D.  Vijayaraghavan 
University  of  Toledo 
Toledo,  Ohio 


July  1995 


Prepared  for 

Lewis  Research  Center 

Under  Cooperative  Agreement  NCC3-291 


U.S.  ARMY 


National  Aeronautics  and 

Space  Administration  RESEARCH  LABORATORY 


An  Efficient  Numerical  Procedure  for 
Thermohydrodynamic  Analysis 
of  Cavitating  Bearings 


D.  Vijayaraghavan* 

Dept  of  Mechanical  Engineering 
University  of  Toledo 
Toledo,  Ohio  43606 


ABSTRACT 

In  this  paper,  an  efficient  and  accurate  numerical  procedure  to  determine  the  thermo¬ 
hydrodynamic  performance  of  cavitating  bearings  is  described.  This  procedure  is  based  on  the 
earlier  development  of  Elrod  for  lubricating  films,  in  which  the  properties  across  the  film 
thickness  are  determined  at  Lobatto  points  and  their  distributions  are  expressed  by  collocated 
polynomials.  The  cavitated  regions  and  their  boundaries  are  rigorously  treated.  Thermal 
boundary  conditions  at  the  surfaces,  including  heat  dissipation  through  the  metal  to  the 
ambient,  are  incorporated.  Numerical  examples  are  presented  comparing  the  predictions  using 
this  procedure  with  earlier  theoretical  predictions  and  experimental  data.  With  a  few  points 
across  the  film  thickness  and  across  the  journal  and  the  bearing  in  the  radial  direction,  the 
temperature  profile  is  very  well  predicted. 


*NASA  Resident  Research  Associate  at  NASA  Lewis  Research  Center. 


INTRODUCTION 


The  increasing  trend  towards  higher  speed,  higher  performance,  but  smaller  size, 
machinery  has  pushed  the  operating  conditions  of  the  bearings  towards  their  'limit  design'. 
Hence,  for  reliable  prediction  of  the  performance  of  such  bearings,  a  model  which  accounts  for 
all  the  operating  conditions  as  closely  as  possible  is  becoming  increasingly  important.  Since, 
the  lubricant  viscosity  strongly  depends  on  temperature,  the  usual  classical  assumptions  of 
constant  viscosity  or  effective  viscosity  become  untenable.  In  many  cases,  the  variation  in 
viscosity  can  cause  a  profound  impact  on  performance  and  it  can  no  longer  be  considered  to 
have  only  secondary  or  perturbation  effect 

The  evidence  of  growing  interest  among  researchers  in  studying  the  thermal  effects  is 
attested  by  the  recent  survey  by  Khonsari  (1987)  and  the  workshop  on  the  thermal  aspects  of 
fluid  film  tribology  (Pinkus,  1990).  The  methodology  of  inclusion  of  thermal  effects  spans 
from  approximations  by  expressing  the  local  temperature  distribution  by  a  single  value,  de¬ 
coupling  of  the  energy  and  Reynolds  equations  and/or  assuming  the  bounding  surfaces  to  be 
adiabatic,  to  a  full  THD  model  analysis  considering  heat  generation  to  dissipation  through  the 
entire  system.  Each  of  these  methods  has  its  own  merits  and  limitations.  Approximations  are 
aimed  at  faster  solutions  while  THD  models  are  used  for  accurate  solutions  at  the  cost  of 
computer  time. 

The  temperature  variation  in  a  lubricant  depends  significantly  on  the  lubricant  flow  rate 
through  the  entire  bearing.  The  film  rupture  and  reformation  locations,  reverse  flow  and 
viscous  dissipation  and  heat  transfer  to  metal  surfaces  in  the  cavitated  region  are  important  in 
determining  the  temperature  profile.  Although  most  of  the  studies  on  thermal  effects  have 
included  the  variations  in  temperature  across  the  film,  generally  the  cavitation  effects  have  been 
only  superficially  considered.  Jakobsson-Floberg  (1957)  -  Olsson  (1965)  have  developed 
boundary  conditions  (JFO  theory),  which  are  based  on  the  conservation  of  mass  principle  and 
the  assumption  that  the  fluid  flows  in  the  cavitated  region  in  the  form  of  striations  extending 
between  the  cavitation  boundaries  and  across  the  film  thickness.  This  model  is  known  to  have 
deficiencies;  but  nevertheless  provides  considerable  improvement  over  the  earlier  more 
approximate  models  (Brewe,  1986). 

More  recently,  there  have  been  many  contributions  in  the  literature  on  thermal  analysis, 
which  include  the  cavitation  effects.  Boncompain  et  al  (1986)  used  Reynolds  boundary 
conditions  for  both  rupture  and  reformation,  and  for  the  film  streamlets  in  the  cavitated  region 
applied  the  energy  equation  of  the  same  form  as  in  full  film,  except  that  the  pressure  and 
transverse  velocity  were  taken  as  zero.  Ott  and  Paradissiadis  (1988)  applied  JFO  boundary 
conditions  to  determine  the  cavitated  region  and  used  the  same  energy  equation  for  both  full 
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film  and  cavitated  region.  However  in  the  energy  equation,  they  used  mean  value  of  fluid 
velocity  components,  averaged  over  the  film  thickness.  Han  and  Paranjpe  (1990, 1994a)  also 
applied  JFO  boundary  conditions  to  determine  the  cavitated  region  using  the  numerical 
algorithm  developed  by  Elrod  (1981).  They  used  a  complete  energy  equation  for  the  entire 
bearing,  with  the  mass  fraction  of  fluid  in  the  cavitated  region  included  in  the  energy  equation. 
Later,  Paranjpe  (1992)  and  Paranjpe  and  Han  (1994b)  extended  the  analysis  for  dynamically 
loaded  bearings  and  non  Newtonian  effects  of  the  lubricant.  Mittwollen  and  Glienicke  (1990) 
considered  mass  conservation  in  their  study  on  the  thermal  effects  on  the  multi-lobed  bearings 
for  turbo  machinery. 

Although  viscosity  variations  across  the  film  thickness  have  decisive  effect  on  the 
bearing  performance,  viscosity  profile  in  that  direction  is  rather  more  important  than  the 
individual  pjoint  wise  viscosity  values.  This  is  because,  the  coefficients  in  the  hydrodynamic 
equation  are  functions  of  integral  of  viscosity  across  the  film.  By  choosing  a  few  points  across 
the  field,  the  location  of  these  points  being  zeros  of  the  Jacobi  polynomial,  the  profile  can  be 
well  predicted.  Elrod  (1991)  showed  that  the  conformity  using  Legendre  polynomials  at 
Lobatto  points,  is  much  better  than  with  polynomials  based  on  equidistant  points.  If  the  order 
of  polynomial  is  increased,  equispaced  interpolation  may  fail  to  converge  while  the  Legendre 
polynomial  interpolation  becomes  progressively  better.  The  Lobatto  integration  considers  the 
end  point  values  also  and  is  a  highly  accurate  scheme.  Obviously,  lesser  the  number  of  points 
needed  across  the  film,  lesser  the  computational  efforts  required.  Also,  when  the  boundary 
conditions  at  film-metal  interfaces  are  in  terms  of  temperature  gradients,  the  unknown 
boundary  temperatures  are  expressed  in  terms  of  all  internal  Lobatto  point  temperatures  using 
the  polynomial  approximation.  Hence,  the  boundary  conditions  will  have  direct  impact  on  all 
the  internal  points,  instead  of  only  the  neighboring  nodes. 

In  this  paper,  a  numerical  procedure  which  is  both  accurate  and  fast  is  described.  This 
is  achieved  by  sampling  the  transverse  temperature,  flow  and  velocity  values  at  few  chosen 
points  and  interlinking  them  by  constructing  orthogonal  polynomials  to  determine  their  profile. 
Elrod  (1991)  introduced  this  basic  approach  for  films.  Later,  Elrod  and  Vijayaraghavan  (1995) 
extended  the  procedure  to  include  cavitation  effects  also.  In  this  procedure,  universal 
hydrodynamic  and  energy  equations  are  developed  for  the  entire  bearing  and  mass  and  energy 
conservation  across  the  reformation  boundary  is  rigorously  implemented.  A  dual  time  step 
procedure  is  used  to  solve  the  THD  model  for  steady  state  results.  Predictions  of  pressure  and 
temperature  profiles  using  this  procedure  are  compared  with  various  analytical  and 
experimental  data  available  in  the  open  literature. 
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NOMENCLATURE 


c 

Radial  clearance 

P  - 

Bulk  modulus 

- 

Specific  heat 

e  - 

Eccentricity  ratio 

D 

. 

Journal  diameter 

<D  - 

Viscous  dissipation 

g 

- 

Switch  function 

K  - 

Thermal  diffusivity 

h 

Film  thickness 

■ 

Viscosity 

ht 

- 

Heat  transfer  coefficient 

0  - 

Partial  film  content 

K 

- 

Thermal  conductivity 

P  - 

Density 

L 

- 

Bearing  length 

t)  - 

Volume  of  cell 

m 

Mass  flux 

0)  - 

Angular  velocity 

N 

- 

No.  of  Internal  Lobatto  points 

^  - 

Ruidity  (l/p.) 

P 

- 

Legendre  polynomial 

V  - 

Angular  co-ordinate 

P 

- 

Pressure 

c  - 

Transformed  co-ordinate  (2z/h) 

Q 

- 

Flow  rate 

q 

- 

Heat  flux 

Subscripts 

r 

- 

Radius 

a  - 

Ambient 

T 

- 

Temperature 

b  - 

Bearing 

t 

- 

Time 

c  - 

Cavitated  region 

U 

- 

Surface  velocity  in  x  direction 

f  - 

Film 

u 

- 

Fluid  velocity  in  x  direction 

g  - 

Groove 

V 

- 

Velocity  vector  (u  +  v  Cy) 

j  - 

Journal 

V 

- 

Total  velocity  vector  (u  -i-  v  ey+  w  e^) 

L  - 

Lower  wall  (^=  -1) 

V 

- 

Fluid  velocity  in  y  direction 

R  - 

Reference 

w 

- 

Fluid  velocity  in  z  direction 

s  - 

Supply 

X 

- 

Circumferential  co-ordinate  axis 

U  - 

Upper  wall  (1=  +1) 

y 

- 

Axial  co-ordinate  axis 

Co-ordinate  axis  across  the  film  thickness 


Underline 

A,B,C,D,E,F,G,R  &  S 


Legendre  coefficients 

Matrices  (see  Numerical  Formulations) 


THEORETICAL  DEVELOPMENT 


The  continuity,  momentum  and  energy  equations  considered  in  this  development  are  of 
the  following  form: 


Continuity : 

V.  v  =  0 

(1) 

a  av 

Momentum : 

1  =  vp 

(2) 

/ar 

^a?T  /avf 

Energy: 

pcpj-^  +  v.vi] 

=  K  —  +  Uhr- 

a^  \az/ 

(3) 

Basic  assumptions  are:  only  viscosity  varies  with  respect  to  temperature;  density  and  other 
thermal  properties  are  constant;  pressure  is  constant  across  the  film  thickness;  Only  the  film 
conduction  and  velocity  gradients  across  the  film  thickness  are  significant.  All  these 
assumptions  are  generally  accepted  for  such  applications,  but  does  not  impose  any  limitation  of 
this  procedure. 

The  fluidity,  which  is  the  reciprocal  of  viscosity,  is  used  as  a  variable  and  expressed 
in  a  series  of  Legendre  polynomials,  as 

^(c)=s&Pk(c)  (4) 

The  Legendre  coefficients  are  obtained  by  collocation  at  Lobatto  points.  Also,  the  fractional 
film  gap,  z,  is  transformed  into  ^  co-ordinates  using  the  relationship  ^  =  z/(h/2),  z  being 
measured  from  the  mid  film  surface.  The  governing  equations  for  the  Legendre  coefficients 
(Elrod,  1991)  are 

V  =  VL  +  Aj^|((p)d(p  +  B£(p^((p)d(p 
where 

Vu-Vl-|b^i 

A= - ^ - 

2^0 

The  hydrodynamic  equation  is  written  as 
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(7) 


where 


12^  +  6(Vu+VL).Vh.6{Vu-VL)V. 
dt 


=  V.^pl^Vp 


^  3^0 

Here,  the  expression  for  is  true,  irrespective  of  number  of  Lobatto  points. 
The  three  dimensional  energy  equation  is  written  as. 


^+V.  VT  + 
dt 


^  ^  dl^  ^  ^ 


where 


Dt 


VdC 


(8) 


Cavitating  Bearing 


In  the  cavitated  region,  the  pressure  is  taken  to  be  constant  at  cavitation  pressure  and  the 
Couette  flow  is  assumed  to  exist  in  the  form  of  striations  occupying  a  portion  of  the  volume. 
The  striations  extend  between  the  cavitation  boundaries  and  bounding  surfaces  (see  Fig.  1).  In 
a  fashion  similar  to  the  development  of  a  single  'universal  equation'  applicable  for  both  full  film 
and  cavitated  regions  for  isothermal  systems  (Elrod  and  Adams,  1974),  the  hydrodynamic 
equation  for  the  present  application  can  be  written  as. 


12  ^  +  6  (Vu  +  Vl)  .  V(9h)  -  6  (Vu  -  Vl)  V. 


where 


0  =  { 


p/Pc  in  full  film  region 

in  cavitated  region 


(9) 


P 

P 


=  Pc  + 


=  e^ 

36 

P(9-1) 


}  for  0  >  1 
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Equation  (9)  has  only  0  as  independent  variable  and  when  solved  determines  0  distribution  for 
the  entire  bearing  (thereby  pressme  distribution)  at  any  instant  of  time. 

According  to  JFO  theory,  the  liquid  film  in  the  cavitated  region  is  confined  to  the 
striations.  Since,  our  interest  is  to  determine  the  film  temperature  variation  in  the  film,  the 
energy  equation  at  the  cavitated  region  is  basically  for  the  striations  only.  Hence,  Eq.  (8)  is 
applicable  for  the  striation  also.  Since  length  to  width  ratio  of  striation  will  be  large,  Elrod  and 
Vijayaraghavan  (1995)  have  developed  the  formulation  to  evaluate  the  Lagrangian  derivative 
(D^/Dt)  for  the  striations,  using  the  short  bearing  theory  which  involves  some  transverse 
pressure  gradient.  It  has  been  shown  that,  within  the  cavitated  region,  independent  of  striation 
width,  we  can  write  the  flow  gradient  term  as. 


Note  that,  at  ^  =  1  this  expression  reduces  to  -2(9h/9t)  and  at  both  upper  and  lower  walls  the 
Lagrangian  derivative  vanishes,  as  it  should.  Also  note  that  there  is  no  transverse  temperature 
variation  within  the  striation  and  since,  the  striations  are  separated  from  one  another,  there  will 
not  be  any  convective  heat  transfer  in  the  axial  direction.  Since,  we  are  only  analyzing  the  film 
in  the  cavitated  region,  there  is  no  need  to  include  the  striation  width  (0)  or  vary  the  thermal 
properties  in  the  cavitated  region.  However,  when  heat  flux  to  the  boundary  surfaces  are 
determined  in  this  region,  corresponding  striation  width  should  be  considered. 

Now,  we  have  one  governing  hydrodynamic  and  one  energy  equation  applicable  to  the 
complete  bearing. 
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Interface  Conditions 


For  cavitating  bearings,  the  film  reformation  interface  has  to  be  treated  more 
jjgQj’ougly^  Figure  1  shows  the  film  reformation  interface  where  the  striated  film  abruptly 
coalesces  to  form  a  full  film.  On  the  cavitated  side  of  the  interface,  the  film  has  a  Couette  flow 
velocity  distribution,  while  on  the  full  film  side,  the  retarding  effect  of  rising  pressure  reduces 
the  film's  convective  capacity.  Across  the  interface,  the  total  mass  and  energy  must  be 
conserved.  Conservation  of  total  mass  flow  is  automatically  satisfied  when  the  hydrodynamic 
equation,  Eq.  (9),  is  solved.  When  there  is  no  reverse  flow  on  the  film  side  of  the  reformation 
interface,  the  Lagrangian  derivative  (D^/Dt)  adequately  describes  the  fluid  motion  across  the 
interface  and  would  satisfy  the  conservation  of  energy.  However,  due  to  adverse  pressure 
gradient  in  the  film  side,  if  a  reverse  flow  occurs  in  the  film  side  near  the  stationary  wall,  then  a 
more  detailed  consideration  of  the  situation  is  warranted.  Elrod  and  Vijayaraghavan  (1995) 
have  developed  a  detailed  procedure  for  determining  flow  pattern  and  temperature  at  the  film 
side  of  the  interface  for  such  cases.  It  is  reasonable  to  assume  that  across  the  interface,  the 
temperature  associated  with  a  fluid  particle  does  not  change  much.  Thus  across  the  interface,  it 
is  hypothesized  that 

T  =  T(Q) 


where 


0(0  ={ 


udC  in  cavitated  region 
ud^  in  full  film  region 


Within  the  reverse  flow  region,  when  Q  <  0  and  u  <  0,  the  temperature  is  determined  by  the 
energy  equation,  from  the  upstream  values.  After  the  velocity  becomes  positive,  but  Q  <  0, 
then  for  the  same  flow,  Q, 

T(Q,u>0)  =  T(Q,u<0)  (11b) 

Beyond  this  region  (Q  >  0  and  u  >  0),  the  temperature  is  determined  from  the  flow  from 
cavitated  side.  For  the  same  flow,  Q, 

[T(Q)]filin  side  ~  [T(Q)]cav.  side  (11c) 

In  other  words,  T(Q)  is  continuous  across  the  reformation  interface,  but  T(C)  may  be 

discontinuous. 

At  the  film  rapture  interface,  since  0=1  at  both  sides  of  the  separation  front,  no  special 
treatment  of  the  interface  is  needed. 
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Boundary  Conditions 


While  the  cavitation  boundary  conditions  are  built  into  the  governing  equations,  the 
thermal  boundary  conditions  for  the  fluid-solid  interfaces  have  to  be  specified.  There  are  a 
variety  of  boundary  conditions  that  can  be  imposed,  based  on  the  total  thermal  system  being 
considered.  When  the  thermal  analysis  within  the  film  is  of  importance,  the  heat  flow  to  the 
boundaries  of  the  film  are  characterized  by  simplified  approximations,  namely,  constant 
temperature,  adiabatic  or  no  net  heat  flow  boundary.  In  a  full  THD  model,  in  addition  to 
considering  the  heat  transfer  within  the  fluid,  heat  dissipation  to  the  ambient  through  the 
boundary  surfaces  must  also  be  determined.  In  this  way,  an  actual  thermal  system  is 
considered  for  the  analysis  and  hence  the  predictions  will  be  more  realistic. 

At  the  fluid-solid  interfaces  the  temperatures  and  the  fluxes  must  be  continuous.  At  the 
film-bearing  interface  these  conditions  yield; 


[TJz  =  h^  -  [Tlr  =  rbi 


&rKf 


at 

dz 


■m 


(12) 


where 

e  for  0  <  1 

In  the  cavitated  region,  the  heat  can  be  transferred  to  the  bearing  surface  only  at  the  striation 
surfaces.  Hence,  the  parameter  6j  is  introduced,  as  0  indicates  the  striation  width. 

The  continuity  of  temperature  and  flux  at  the  fluid-journal  interface  is  also  imposed  in 
similar  manner,  except  that,  since  the  journal  is  in  motion,  the  cyclic  variation  of  journal 
surface  temperature  can  be  assumed  to  be  constant  in  the  circumferential  direction  (Dowson  et 
al,  1966).  Therefore,  the  heat  flux  to  the  journal  surface  can  be  averaged  over  the 
circumference.  The  net  heat  flux  is  determined  by  integrating  the  film  temperature  gradient  in 
the  following  manner  (Khonsari  and  Beaman,  1985). 


[T]z  =  -h/2  =  [T]r  =  rj 


Kf 
27C  r 


(13) 


The  outer  surfaces  and  the  lateral  surfaces  of  the  bearing  transfer  heat  to  the  ambient  through 
free  convection,  depending  upon  the  prevailing  heat  transfer  coefficient 
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In  the  case  of  grooved  journal  bearings,  the  hot  oil  earned  over  mixes  with  the 
relatively  cooler  supply  lubricant.  Assuming  instantaneous  mixing  at  the  groove,  the 
temperature  of  the  lubricant  leaving  the  groove  is  obtained  using  a  simple  heat  balance.  Often, 
the  instantaneous  mixing  assumption  can  be  questionable  in  real  situations.  Heshmat  and 
Pinkus  (1986)  and  Mitsui  et  al  (1983),  among  others,  have  tested  a  number  of  different 
bearings  and  have  developed  empirical  equations  to  estimate  the  mixing  coefficient  used  to 
determine  the  mixed  lubricant  temperature.  However,  for  the  numerical  examples  presented 
here,  only  simple  heat  balance  is  used  to  determine  the  oil  mixture  temperature  at  the  groove. 


Heat  Conduction  Equations 


The  heat  transfer  within  the  circular  bearing  shell  is  governed  by  the  following  heat  conduction 
equation,  written  in  cylindrical  co-ordinate  system. 


[^Tb  ,  1  axi,  1  S^Tt 


(14) 


In  the  case  of  a  rotating  journal,  the  journal  temperature  is  considered  to  be  independent  of  the 
X  direction  and  the  relevant  term  drops  out. 

Since,  the  hydrodynamic  and  energy  equations  are  set  up  in  unsteady  form,  the  above 
heat  conduction  equations  are  also  written  in  unsteady  form.  Although  in  the  present  paper, 
only  steady  state  cases  are  considered,  the  development  of  this  procedure  is  aimed  at  analyzing 
dynamically  loaded  applications.  Also,  for  steady  state  cases,  the  time  derivative  term  in  these 
governing  equations  provides  time  dependent  numerical  dissipation  and  has  a  stabilizing  effect 
on  the  iteration  process.  Here,  the  time  step  is  used  as  a  control  parameter  for  the  rate  of 
convergence  towards  asymptotic  steady  state  conditions 

Equations  (9),  (8),  (10)  and  (14),  along  with  the  appropriate  boundary  conditions,  are 
solved  simultaneously  to  determine  the  temperature  distribution  for  the  whole  system. 


10 


NUMERICAL  FORMULATIONS 
Matrices 


Elrod  (1991)  has  indicated  the  procedure  to  set  up  the  required  matrices.  For  the  sake 
of  continuity,  the  overall  approach  is  mentioned  here.  It  is  convenient  to  set  up  the  fluidity  in  a 
vector  form  and  determine  the  velocity,  flow  and  temperature  using  matrices.  Using  matrix 


notation  we  can  wnte 

N+l  N+1 

^(Ck)  =  I  ^Pn(Ck)  =  I  Ckn^ 

0  0 


(15) 


where  is  the  Legendre  polynomial  of  order  n,  ^  is  the  Legendre  coefficient,  is  the  ^ 
value  at  k*  Lobatto  point  and  N  is  the  number  of  internal  Lobatto  points.  Therefore, 

^  and  ^  (16) 

Similarly,  we  can  also  write, 

T  =  Cl  andT  =  C*^T  (17) 

The  velocity  profile  according  to  Eq.(5)  can  be  determined  from, 

V  =  Vl  +  AF^  +  B  G4  (18) 

The  matrices  F  and  G  can  be  obtained  by  evaluating  the  integrals  in  Eq.(5)  at  Lobatto  points. 
Similarly,  flow  rate  obtained  by  integrating  the  velocity  profile,  is  computed  using 

£vdC  =  Vl(i  +  C)  +  AR^  +  B 

The  temperature  gradients  in  the  C  direction  is  obtained  by  successively  differentiating  Eq.(17), 
and  can  also  be  written  in  matrix  form  as, 

—  =D1  and  ^  =  ET 

(20) 

Once  the  number  of  internal  points  is  known,  the  Lobatto  point  locations,  weighting 
factors  and  the  Legendre  polynomials  can  be  determined  using  standard  mathematical 
handbook;  for  example  refer  Abramovitz  and  Stegun  (1965).  Then,  the  matrices  C,  C  \  D,  E, 
F,  G,  R  and  S  can  be  easily  determined.  Note  that,  these  matrices  need  to  be  computed  only 
once. 
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Cavitation  Algorithm 


The  'universal'  hydrodynamic  equation,  Eq.(9),  is  written  in  finite-difference 
formulation  using  the  Elrod  cavitation  algorithm  (Elrod,  1981),  appropriately  modifying  the 
terms  for  the  present  form  of  the  equation.  The  cavitation  algorithm  is  solved  using  the  two 
step  alternate  direction  implicit  (ADI)  method,  in  which  the  time  step  is  split  into  two  halves 
and  for  each  half,  0  values  in  one  co-ordinate  direction  is  implicitly  solved.  For  each  sweep, 
this  results  in  only  tridiagonal  system  of  equations  to  be  simultaneously  solved.  The  tridiagonal 
matrices  can  be  efficiently  inverted  using  Thomas  algorithm  (refer  Anderson  et  al,  1984). 
When  the  process  is  completed,  6  value  distribution  is  updated  for  next  time  step. 

Energy  Equation 

The  energy  equation,  Eq.(8),  is  implicitly  solved  using  three  dimensional  ADI  method, 
developed  by  Douglas  and  Gunn  (1964)  (also  refer  Anderson  et  al,  1984).  This  procedure, 
which  is  an  extension  of  Crank-Nicholson  scheme  to  three  dimensional  form,  is  second-order 
accurate  in  time  and  space,  and  is  stable.  The  solution  for  one  time  interval  is  obtained  in  three 
steps.  At  each  step,  one  direction  is  treated  implicitly  while  the  derivatives  in  other  directions 
are  evaluated  for  known  values  of  temperature.  With  the  present  formulation  of  the  energy 
equation,  during  the  sweep  in  C  direction,  (N  x  N)  matrix  has  to  be  inverted,  N  being  the 
number  of  internal  Lobatto  points;  for  the  other  two  directions  only  tridiagonal  matrices  have  to 
be  inverted. 

For  the  full  THD  model,  the  temperature  distribution  in  the  fluid,  bearing  and  journal  is 
to  be  determined.  In  order  to  efficiently  and  simultaneously  determine  the  temperature 
distribution  in  the  bearing/runner  metals  and  also  to  effectively  couple  the  boundary  conditions 
at  fluid-metal  interfaces,  Lobatto  points  are  positioned  within  both  the  bearing  and  journal 
domains.  Now,  the  bearing  and  journal  temperature  gradients  in  C  direction  can  be  expressed  in 
matrix  form,  similar  to  Eq.(20).  Proper  care  must  be  taken  to  index  all  the  Lobatto  points, 
since,  different  equations,  boundary  conditions  and  variables  are  being  used  in  each  domain. 

When  there  are  regions  of  reverse  flow  in  the  bearing,  computational  instabilities  could 
be  encountered.  The  negative  velocity  hinders  the  solution  of  energy  equation,  which  is  being 
treated  as  initial  value  problem  and  marched  forward.  Also,  the  convection  terms  in  the  energy 
equation  are  of  first  order.  Central  differencing  of  these  terms  sometimes  causes  numerical 
instability  due  to  lack  of  diagonal  terms.  To  alleviate  both  problems  and  to  preserve  numerical 
stability,  one  sided  upwind  differencing  is  employed  for  the  convection  terms.  For  example. 
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(21) 


u  =  abs  (ui)  igx  =  sign  (ui) 

ox|i  Ax 

Similarly  for  convection  in  transverse  (y)  direction. 

Boundary  Conditions 

When  the  boundary  smfaces  are  held  at  constant  temperature,  only  the  temperature 
distribution  at  internal  Lobatto  points  need  be  determined.  However,  when  one  or  more 
boundaries  have  other  boundary  conditions  involving  temperature  gradients,  those  boundary 
temperatures  are  also  unknown.  Because  of  our  polynomial  approximation  procedure,  it  is 
possible  to  express  the  boundary  temperatures  and  thus  the  boundary  conditions  in  terms  of  all 
the  internal  Lobatto  point  temperatures.  In  the  THD  model  temperature  at  all  the  four 
boundaries,  two  fluid-metal  interfaces,  metal-ambient  interface  and  journal  center  or  inside 
diameter,  are  unknown.  It  is  still  possible  to  express  these  boundary  conditions  in  terms  of 
internal  Lobatto  point  temperatures  in  the  form  of  simultaneous  equations.  The  temperamre 
distribution  in  the  whole  system  can  be  solved  in  implicit  manner  and  from  the  internal 
temperatures,  boundary  temperatures  can  be  determined.  Though,  this  process  is  somewhat 
tedious,  it  is  relatively  straight  forward.  For  more  details,  interested  readers  are  referred  to 
Villadsen  and  Michelsen  (1978). 

Dual  Time  Step 

The  thermal  response  of  a  fluid  film  is  of  the  order  of  milliseconds,  whereas,  the 
temperature  of  metal  parts  takes  seconds  or  minutes  to  stabilize  due  to  their  large  thermal 
capacity.  Also,  with  the  numerical  procedure  point  of  view,  the  grid  spacing  across  the  film 
thickness  is  of  the  order  of  microns,  while  in  the  metal  the  grid  points  are  in  terms  of 
millimeters  apart.  When  the  whole  system  is  solved  simultaneously,  these  differences  can 
cause  numerical  problems.  Khonsari  and  Beaman  (1985)  reported  that  when  the  intermediate 
results  of  the  metal  conduction  equation  are  directly  used  as  boundary  conditions  for  energy 
equation,  the  solution  becomes  unstable.  They  used  all  equations  in  steady  state  form  and  the 
Reynolds  and  energy  equations  were  solved  several  times  before  moving  to  heat  conduction 
equation.  Paranjpe  and  Han  (1994b)  also  discussed  about  the  time  scales  for  thermal  transients 
and  determined  that  the  time  scale  for  thermal  transients  for  journal  and  bearing  to  be  3  to  4 
orders  of  magnitude  greater  than  those  for  fluid. 

To  tackle  this  problem,  a  dual  time  step  procedure  is  utilized  here.  Smaller  time  steps 
(of  the  order  of  lO'^-lO"^  sec)  were  used  to  solve  fluid  film  equations  and  larger  time  steps  (of 
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the  order  of  10^-10'^  sec)  were  used  in  bearing  and  runner  metal  domains,  and  all  the 
equations  are  solved  simultaneously.  For  steady  state  problems,  the  system  eventually  attains 
an  equilibrium  conditions  with  the  surroundings  and  this  dual  time  scale  should  not  matter. 
With  this  dual  time  march,  we  are  only  accelerating  the  temperature  changes  in  metal 
compatible  with  the  temperature  changes  in  fluid.  This  procedure  is  found  to  be  very  stable 
and  relatively  faster  in  convergence.  The  final  solution  was  also  found  to  be  independent  of  the 
time  step.  When  the  time  step  used  for  conduction  through  the  metal  is  not  compatible  with  the 
time  step  used  in  fluid  domain  and  the  solution  diverges,  then  the  time  steps  are  adjusted  and 
the  analysis  is  repeated. 

Solution  Procedure 

When  the  number  of  Lobatto  points  is  known,  their  locations  and  all  the  matrices  are 
readily  computed.  Initially,  the  whole  bearing  clearance  is  assumed  to  be  filled  with  lubricant 
and  the  internal  temperatures  are  set  as  a  function  of  boundary  temperatures.  Fluidity  values  at 
all  points  are  determined  according  to  the  latest  available  temperatures  and  the  viscosity- 
temperature  relationship.  The  hydrodynamic  equation  is  solved  using  the  cavitation  algorithm 
to  determine  the  0  distribution  and  thereby  the  pressure  and  mass  flow  rate.  Using  the  flow 
vectors  A  and  B  and  other  matrices,  three  dimensional  velocities,  flows  and  viscous  dissipation 
at  every  Lobatto  point  are  determined.  Across  the  film  reformation  front,  the  interface 
conditions  as  described  above  are  enforced.  The  energy  equation  is  implicitly  solved.  After 
each  step,  the  boundary  temperatures  are  updated.  After  each  time  step  the  convergence  in  6 
and  T  are  checked  and  if  not  adequate,  the  process  is  repeated  by  updating  the  fluidity  values 
and  switch  function,  g,  values.  In  the  case  of  grooved  journal  bearings,  at  every  time  step  the 
mixed  lubricant  temperatm^  is  determined  using  simple  heat  balance.  The  nodes  in  the  groove 
location  are  assumed  to  be  at  constant  supply  pressure  and  at  constant  mixed  lubricant 
temperature. 

The  numerical  code  is  written  for  finite  bearings  and  is  vectorized  for  use  on  Cray 
XMP/YMP  computers.  For  the  numerical  examples  discussed  in  this  paper,  convergence 
accuracies  of  10'^  for  0  and  10'^  for  temperature  are  imposed.  However,  it  is  found  that  the 
convergence  criteria  for  0  is  more  stringent  and  when  the  steady  state  conditions  are  achieved, 
the  temperature  convergence  is  below  10  .  Various  grid  sizes  and  Lobatto  points  were  tried 
and  it  was  found  that  above  41  nodes  in  the  circumferential  direction,  depending  upon  L/D,  7 
to  1 1  nodes  along  half  the  bearing  length  in  axial  direction  and  3  Lobatto  points,  the  results  are 
very  similar;  but  the  computer  time  requirement  was  becoming  larger.  In  order  to  make  use  of 
vectorization  of  Cray  computers,  49  nodes  in  the  circumferential  direction,  9-11  nodes  axially. 
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3  internal  Lobatto  points  across  film  and  journal  and  5  internal  Lobatto  points  across  bearing 
radius  are  chosen.  The  computer  time  requirement  is  primarily  dependent  on  the  type  of 
problem  being  solved.  For  a  typical  finite  journal  bearing,  to  solve  a  direct  problem  with 
eccentricity  specified,  with  adiabatic  boundary  conditions,  the  CPU  time  on  a  Cray  XMP 
computer,  for  the  above  grid  arrangement,  was  15-20  seconds  and  for  a  typical  THD  model 
including  journal  and  bearing,  the  CPU  time  requirement  was  100-200  seconds. 

NUMERICAL  EXAMPLES 

With  the  present  procedure,  several  cases  are  analyzed  and  the  results  are  compared 
with  the  earlier  analytical  predictions  and  experimental  data  available  in  literature.  The  bearing 
data  are  shown  in  Table  1. 

Ott  and  Paradissiadis  (1988)  and  later  Han  and  Paranjpe  (1990)  analyzed  the 
temperature  distribution  in  an  axial  grooved  journal  bearing  wherein  the  journal  and  oil  inlet 
groove  are  maintained  at  constant  oil  inlet  temperature  and  the  bearing  surface  is  adiabatic.  The 
pressure  and  temperature  contours  at  the  axial  symmetry  plane  predicted  with  just  3  internal 
Lobatto  points  are  compared  with  the  predictions  of  Han  and  Paranjpe  (1990),  as  shown  in 
Figs.  2  and  3,  respectively.  The  trends  in  these  figures  are  similar  to  the  earlier  predictions. 
Both  the  present  prediction  and  those  of  Han  and  Paranjpe  have  the  same  pressure  profile  with 
a  peak  pressure  of  18  MPa,  compared  to  the  value  of  14.5  MPa  predicted  by  Ott  and 
Paradissiadis.  The  temperature  contours  have  some  distinct  differences.  Both  Ott  and 
Paradissiadis  and  Han  and  Paranjpe  predicted  the  fluid  temperatures  in  the  cavitated  region  to 
be  fairly  constant  at  the  same  C  locations,  with  the  maximum  value  being  100*^C  and  120®C 
respectively.  With  the  present  procedure,  as  can  be  seen  from  Fig.  (3),  fluid  temperatures  in  the 
cavitated  region  increase  at  the  same  ^  locations,  with  the  maximum  value  being  130®C  and 
fluid  with  lower  temperature  is  closer  to  the  constant  temperature  moving  journal.  In  the 
cavitated  region,  the  film  thickness  inaeases  and  the  fluid  velocity  component  in  the  ^  direction 
will  have  to  be  significant  in  order  to  fill  in  the  gap,  towards  the  hotter  adiabatic  boundary 
surface.  Therefore,  it  is  reasonable  to  expect  the  temperature  to  rise  in  this  region,  as  predicted 
by  the  present  procedure.  The  load  capacity,  as  per  the  present  prediction  is  328.4  kN  as  against 
273.2  kN  obtained  by  Ott  and  Paradissiadis. 

Ferron  et  al  (1983)  conducted  experiments  to  measure  the  temperature  variations  in  a 
journal  bearing  and  compared  the  data  with  their  theoretical  predictions.  They  determined  that 
due  to  differential  dilatation  of  the  steel  journal  and  bronze  bearing,  the  nominal  radial  clearance 
of  145  mm  at  20°C  became  152  mm  at  45°C.  Later  Boncompain  et  al  (1986)  in  their  analysis 
of  the  same  bearing,  in  addition,  approximated  the  distorted  bearing  internal  surface  to  be  a 
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circle  and  determined  that  the  nominal  radial  clearance  of  145  mm  and  eccentricity  ratio  of  0.8 
at  20^C  is  equivalent  to  157  mm  and  0.74  respectively  at  45  C.  Khonsari  and  Wang  (1991) 
also  performed  analysis  on  this  bearing,  by  including  thermal  expansion  for  both  journal  and 
bearing  and  elastic  deformation  of  bearing  and  determined  that  the  collective  effect  of  these 
factors  yielded  satisfactory  comparisons  with  the  experimental  data.  Figure  4  is  the  comparison 
of  temperature  distribution  at  the  bearing  symmetry  plane,  predicted  using  3  internal  Lobatto 
points  each  across  film  and  journal  radius  and  5  internal  Lobatto  points  across  bearing 
thickness  with  the  predictions  of  Boncompain  et  al  (1986).  Obviously,  the  plots  are  not  to 
scale,  with  the  dimensions  of  the  film  region  highly  exaggerated,  so  that  temperature  variations 
within  the  film  can  be  clearly  seen.  The  general  trend  in  the  temperature  contours  predicted  by 
both  procedures  is  similar.  Both  analyses  predict  a  peak  temperature  of  about  52  C.  However, 
the  present  analysis  predicts  the  maximum  fluid  temperature  around  the  film  rupture  region, 
whereas,  Boncompain  et  al  predict  this  occurrence  in  the  cavitated  region.  In  the  THD  model, 
since  heat  is  being  dissipated  to  ambient  through  the  bearing  metal,  the  bearing  inside  surface  is 
cooler  than  the  fluid  temperature  in  this  region.  Hence,  the  maximum  temperature  should  occur 
in  the  vicinity  of  minimum  film  thickness,  where  the  viscous  dissipation  is  the  largest.  Notice 
that,  the  oil  film  near  the  supply  groove  is  cooler  than  the  bearing/joumal  metal  temperature.  In 
this  case,  much  of  the  journal  metal  is  at  about  46°C.  The  isotherms  in  the  bearing  metal  are 
primarily  in  the  radial  direction. 

We  now  compare  predicted  bearing  temperatures  with  values  determined  in  various 
experiments.  Figures  5  and  6  compare  the  present  predictions  of  pressure  and  temperature 
profile  around  the  circumference  of  the  bearing  at  the  symmetry  plane  with  the  experimental 
data  of  Ferron  et  al  (1983)  when  the  bearing  is  operating  at  2(XX)  rpm  under  4  kN  load  and  at 
4000  rpm  under  6  kN  load.  The  predicted  pressure  profile,  considering  only  the  differential 
dilatation,  agrees  very  well  with  the  data  and  the  predicted  bearing  inner  surface  temperature 
variation  is  close  to  the  experimental  data,  particularly  in  the  cavitated  region.  Since,  the  active 
part  of  the  thermocouples  are  located  near  the  film-bearing  interface,  the  fluid-metal  interface 
temperatures  are  used  for  this  comparison.  As  can  be  seen  from  Fig.  4,  the  temperature  profile 
in  the  bearing  is  predominantly  in  the  radial  direction  and  temperature  difference  between 
inside  and  outside  surfaces  are  only  about  2-4*^C.  Nondimensional  side  leakage  flow  rate 
(Q/LcU)  predicted  by  the  present  procedure  are  compared  with  the  experimental  results  of 
Ferron  et  al  (1983)  for  a  wide  range  of  eccentricity  ratios  and  at  two  different  speeds  (2(X)0 
rpm  and  4000  rpm)  in  Fig.  7.  The  agreement  is  excellent.  Also  note  that,  for  an  isothermal 
analysis,  irrespective  of  the  speed,  the  nondimensional  flow  rate  will  be  the  same.  The  thermal 
effects  on  the  flow  rate  in  the  bearing  are  clearly  brought  out  in  this  comparison.  Figure  7  also 
provides  plots  of  power  loss  at  these  conditions.  There  is  a  significant  increase  in  the  power 
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loss  due  to  increased  speed,  but  it  does  not  vary  much  due  to  variation  in  eccentricity  ratio  or 
applied  load. 

Mitsui  (1987)  conducted  a  number  of  experiments  on  circular  journal  bearings.  The 
bearing  with  an  L/D  value  of  0.7,  had  one  axial  groove  of  10®  arc,  extending  85%  of  length, 
located  at  the  crown  of  top  pad.  Figure  8  compares  the  predictions  of  bearing  outside 
temperature  with  the  experiments  conducted  at  3.73  kN  and  1750  rpm,  with  two  different 
lubricants,  transformer  oil  and  #140  turbine  oil.  The  predicted  temperature  variation  is  close  to 
the  experimental  data,  particularly  in  the  cavitated  region.  The  higher  the  viscosity  of  lubricant 
used,  the  larger  the  temperature  rise.  Though  the  viscosity  of  turbine  oil  is  more  than  six  times 
that  of  transformer  oil  at  supply  temperature,  the  effective  viscosity  in  the  film  is  not  that  much 
different,  due  to  the  larger  temperature  rise. 

Lund  and  Tonnesen  (1984)  conducted  experiments  on  a  two  axial  groove  journal 
bearing,  each  groove  extending  10®  circumferentially  and  half  the  length  axially.  Figure  9 
compares  the  circumferential  profile  of  the  bearing  wall  temperature  at  the  bearing  symmetry 
plane,  at  the  same  applied  load  was  5.6  kN  and  when  the  journal  speed  3500  rpm  and  5000 
rpm.  The  predicted  temperature  variation  is  in  close  agreement  with  the  experimental  data, 
except  upstream  of  the  second  groove.  Fitzgerald  and  Neal  (1992)  also  conducted  experiments 
on  two  axial  groove  journal  bearings,  each  groove  extending  28.7®  circumferentially  and  80% 
of  the  bearing  length.  Figure  10  provides  a  comparison  of  the  circumferential  profile  of  the 
bearing  wall  temperature  at  the  bearing  symmetry  plane,  when  the  applied  load  was  9.43  kN  at 
the  journal  speed  of  8000  rpm;  but  for  different  UD  ratios  of  0.5  and  1 .0.  Again,  the  predicted 
temperature  variation  is  in  close  agreement  with  the  experimental  data,  except  upstream  of  the 
second  groove.  In  the  twin  axial  groove  bearings,  although  there  is  no  significant  supply  of 
lubricant  through  the  second  groove,  the  groove  is  assumed  to  be  flooded.  Hence,  at  the 
upstream  of  the  second  groove  the  striation  width  is  comparatively  larger.  Since,  the  heat 
transferred  to  bearing  metal  in  the  cavitated  region  is  at  the  striation  surface,  more  heat  is 
dissipated  to  the  ambient  at  this  region.  Therefore,  the  predicted  metal  temperature  in  this 
region  rises  a  bit  before  dropping.  In  the  experiment,  there  could  be  mixing  of  hot  oil  with  the 
cooler  lubricant  causing  the  fluid  temperature  to  reduce  at  the  second  groove,  as  can  be  seen  in 
Fig.  9.  Considering  this,  the  rise  in  temperature  of  the  fluid  at  the  upstream  of  second  groove 
is  compatible  with  the  predictions.  Such  a  situation  is,  of  course,  difficult  to  treat  analytically. 
Knight  and  Niewiarowski  (1990)  also  predicted  the  temperature  profile  for  the  above  cases 
using  a  gas  bubble  model  and  obtained  good  agreement  with  the  experimental  data  in  the 
cavitated  region.  The  gas  bubble  model  assumes  that  the  film  is  continuous  in  the  cavitated 
region,  but  the  effective  density  and  viscosity  of  the  film  are  varied  according  to  the  gas 
fi^ction.  In  the  present  approach,  the  film  occupies  only  a  part  of  the  volume  in  the  cavitated 
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region,  corresponding  to  the  mass  fraction  and  only  the  film  is  considered  for  thermal  analysis. 
Ma  and  Taylor  (1994)  developed  a  seperation  cavitation  model  to  account  for  the  oil  back  flow 
from  the  groove  into  the  cavitated  region  and  were  able  to  determine  the  temperature  fade. 
This  present  model  also,  by  conserving  mass  flow,  does  take  into  account  back  flow  in  the 
groove  and  correctly  determines  the  temperature  fade  in  the  cavitated  region  at  the  trailing 
groove. 

In  general,  the  present  predictions  compare  very  well  with  various  analytical  and 
experimental  data.  While  comparing  the  predicted  results  with  experimental  data,  it  should  be 
recognized  that  there  will  be  number  of  variables  in  the  experiment  that  can  not  be  easily 
simulated  theoretically.  To  name  a  few:  operating  clearance,  film  shape,  mixing  at  the  groove, 
heat  transfer  coefficient,  heat  dissipation  to  other  components  in  the  test  rig  etc.  Considering  all 
these  factors,  it  can  be  said  that  the  present  procedure  predicts  the  temperature  profiles  rather 
well. 

CONCLUSIONS 

An  efficient  numerical  procedure  to  determine  the  thermohydrodynamic  performance 
of  the  cavitating  bearings  is  developed  considering  mass  and  energy  conservations  and  more 
accurately  treating  the  cavitated  region  and  its  boundaries.  This  procedure  requires  few  points 
across  the  film  thickness  strategically  located  and  the  profile  is  expressed  using  collocated 
polynomials.  The  details  of  the  theoretical  and  numerical  formulations  are  presented.  The 
predictions  using  this  procedure  are  compared  with  various  numerical  and  experimental  data, 
and  lead  us  to  the  following  major  conclusions: 

(i)  The  predictions  are  in  good  agreement  with  the  experimental  data  and  predictions  by 
others. 

(ii)  With  just  3  internal  points  across  the  film  (and  across  the  bearing  shell  and  journal 
radius),  the  temperature  profile  can  be  well  predicted. 

(iii)  Although  some  matrix  operations  are  required  with  this  procedure,  the  computer  time 
requirement  is  very  small. 

(iv)  The  approach  gives  a  logical  treatment  of  mass  and  energy  conservation  across  the 
cavitation  boundaries  and  determines  three  dimensional  viscosity,  velocity,  flow,  viscous 
dissipation  and  temperature  profiles. 

(v)  With  a  dual  time  step  procedure,  the  temperature  distributions  across  the  film  and  the 
metals  are  simultaneously  determined,  resulting  in  faster  steady  state  solutions. 
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Figure  1  .—Schematic  diagram  of  a  cavitating  bearing. 
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Figure  2.— Comparison  of  pressure  contours  at  symmetry 
plane,  (a)  Present,  (b)  Han  and  Paranjpe,  1990. 
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Figure  3. — Comparison  of  temperature  contours  at  symmetry 
plane,  (a)  Present,  (b)  Han  and  Paranjpe,  1990. 
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Figure  5. — Comparison  of  pressure  variation  at  symmetry 
plane  (Perron  et  al.,  1983). 


Present  Expt 


6  kN,  4000  rpm  -  O 

4  kN,  2000  rpm - A 


Angle  from  groove  CL,  deg 


Figure  6. — Comparison  of  bearing  inside  surface  temper¬ 
ature  variation  at  symmetry  plane  (Perron  et  al.,  1983). 
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Figure  7. — Comparison  of  flow  rate  and  power  loss 
(Perron  et  al.,  1 983). 
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Figure  8. — Comparison  of  bearing  outside  temperature  at 
symmetry  plane  (Mitsui,  1987). 


5000  - 

rpm  - 

3500  Load:  5.6  kN 
rpm 


A  A 


60 

-V  A  F'I'TI 

rupture^  ^Groove 


o> 


A  \  \ 

A  aY 


0  90  1  80  270  361 

Angle  from  groove,  deg 

Figure  9. — Comparison  of  bearing  outside  temperature  at 
symmetry  plane  (Lund  and  Tonnesen,  1984). 
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Figure  1 0. — Comparison  of  bearing  outside  temperature  at 
symmetry  plane  (Fitzgerald  and  Neal,  1992). 
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